from pylab import *
import plotly.graph_objects as go
import functools
import itertools4 Cálculo de CCM para una cadena quiral C3
En este capítulo se estudiará la quiralidad de los modos normales de una cadena quiral con simetría \(C_{3}\). Los criterios usados son la polarización de fonón, el CCM y el pseudoescalar de momento angular.
Debajo, se definen varias constantes a emplear. Dentro de estas están las masas de los sitios, la longitud de los lados del ángulo vista desde arriba, las constantes de resorte y el ángulo de los enlaces respecto al plano x - y. Además, se definen las posiciones los tres átomos dentro de la celda unitaria, y dos átomos de celdas vecinas. También, se defienen dos rotaciones que serán usadas de forma constante.
# Masas de los distintos sitios.
m1 = 1
m2 = 1
m3 = 1
M = [m1,m2,m3]
# Matriz de resortes a lo largo de x.
KL = 1
KT1 = 0.05
KT2 = 0.25
# Tensor de constante de red en configuración usual.
Kx = array([[KL,0,0],
[0,KT1,0],
[0,0,KT2]],dtype=complex)
# Longitud de los enlaces proyectados al plano
a = 1
# Ángulo de los enlaces respecto al plano x - y.
φxy = pi/3
# Longitud de la celda unitaria.
c = 3 * a * tan(φxy)
avec = array([0,0,c])
#Rotación de φ respecto al plano x - y.
def T(φ):
return array([[cos(φ) ,0. ,sin(φ)],
[0. ,1. ,0],
[-sin(φ),0. ,cos(φ)]])
#Rotación respecto al eje z.
def U(θ):
return array([[cos(θ),-sin(θ),0],
[sin(θ),cos(θ),0],
[ 0, 0,1]] )
# Posiciones de los átomos en la celda unitaria de la cadena.
r1 = array([1,0,0])
r2 = dot(U(2*pi/3),r1) + array([0,0,c/3])
r3 = dot(U(4*pi/3),r1)+ array([0,0,2*c/3])
# Posición de los átomos en la celda vecina.
r1plus = r1 + avec
r3minus = r3 - avec
positions = array([r1,r2,r3])
positionsnn = array([r3minus,r1plus])
X,Y,Z = positions.T
Xnn,Ynn,Znn = positionsnn.TDebajo, se muestra una imagen de la cadena.
DATA = [go.Scatter3d(x=X, y=Y, z=Z,mode='markers',marker_color = "blue",showlegend=False)]
DATA.append(go.Scatter3d(x=Xnn, y=Ynn, z=Znn,mode='markers',marker_color = "red",showlegend=False))
DATA.append(go.Scatter3d(x= [X[0],Xnn[0]], y=[Y[0],Ynn[0]], z=[Z[0],Znn[0]],mode='lines',line=dict(color='red'),showlegend=False))
DATA.append(go.Scatter3d(x= [X[0],X[1]], y=[Y[0],Y[1]], z=[Z[0],Z[1]],mode='lines',line=dict(color='blue'),showlegend=False))
DATA.append(go.Scatter3d(x= [X[2],X[1]], y=[Y[2],Y[1]], z=[Z[2],Z[1]],mode='lines',line=dict(color='blue'),showlegend=False))
DATA.append(go.Scatter3d(x= [X[2],Xnn[1]], y=[Y[2],Ynn[1]], z=[Z[2],Znn[1]],mode='lines',line=dict(color='red'),showlegend=False))fig = go.Figure(data=DATA)
fig.show()4.1 Cálculo de CCM de la estructura
Para calcular el CCM de la estructura se sigue la definición del capítulo pasado. Sin embargo, para hacer este cálculo sólo se consideran reflexiones por planos que contengan al eje z. Esto pues sólo nos interesa la quiralidad respecto a este eje. Aún así, se define una rotación respecto al plano x - y que será utilizada más adelante en el cálculo de la matriz dinámica.
r2 - dot(U(2*pi/3),r1)array([0. , 0. , 1.73205081])También se definen la matriz identidad en 3D y una reflexión respecto al plano xz. Además, se define una funciín que transforma operadores. Esta función se usará para rotar el plano respecto al cual se realiza la rotación.
#Matriz identidad.
Id = array([[1,0,0],
[0,1,0],
[0,0,1]] )
#Reflexión respecto al plano x - z.
σy = array([[-1,0,0],
[0,1,0],
[0,0,1]] )
def TensT(O,A):
"""Función que aplica transforma por un operador O a un tensor A."""
return dot( inv(O), dot(A,O) )Ya con todas las funciones a usar bien definidas, se hace el cálculo de la CCM De la estructura. Para ello, primero se calculan las coordenadas respecto al centro de masa.
# Centro de masa
RCM = (m1*r1 + m2*r2 + m3*r3)/(m1 + m2 + m3)
# Coordenada Relativa 1
qrcm1 = sqrt(m1)*(r1 - RCM)
# Coordenada Relativa 2
qrcm2 = sqrt(m2)*(r2 - RCM)
# Coordenada Relativa 3
qrcm3 = sqrt(m3)*(r3 - RCM)
# Lista con las coordenadas relativas.
Qrcm = [qrcm1,qrcm2,qrcm3]Queremos calcular la CCM tal que el plano de inversión elegido maximize la superposición de la estructura con la estructura invertida. Para ello, se define una función que calcula el CCM de la estructura tomando una reflexión respecto a un plano de entrada.
def CCMS(σ):
num = 0
den = 0
for q in Qrcm:
num = num + dot(q,dot(Id + σ,q))
den = den + dot(q,q)
return 1 - (num/(2*den))Ya hecho esto, se calcula el CCM de la estructura. Notamos que para esta cadena el CCM es igual sin importar que plano de reflexión se use. Esto se ve en la gráfica mostrada debajo.
θ = linspace(0,pi,300)
σ = []
for ang in θ:
σ.append(TensT(U(ang),σy))
σ = array(σ)
CCM_estructura = []
for ref in σ:
CCM_estructura.append(CCMS(ref))
print("La CCM de la estructura es", str(np.min(CCM_estructura)))La CCM de la estructura es 0.16666666666666663Labels = [r"0",r"π/2",r"π"]
Ticks = concatenate([[0],[pi/2],[pi]])
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,4.5))
ax.plot(θ,CCM_estructura)
ax.set_xticks(Ticks)
ax.set_ylim([0,1])
ax.set_xticklabels(Labels,fontsize = 20);
ax.set_ylabel("CCM",fontsize = 20)
plt.show()
4.2 Cálculos del CCM de los modos vibracionales
Ahora calculamos el CCM asociado a los modos normales. Para ello, debemos también construir un programa que obtenga las frecuencias de dichos modos, así como una función que de el CCM de un modo dado.
Primero, se define una función que regresa el valor máximo del producto \(\langle Q| Q \rangle\) para cierto plano de reflexión. Además, se define una función que regresa el valor esperado y que será usada más adelante.
def expectedval(vec,S):
return dot(np.conj(vec).T,dot(S,vec)).real
θ = linspace(0,pi,300)
def CCMMax(q):
max = 0
for ang in θ:
σ = TensT(U(ang),σy)
if abs(expectedval(q,Id + σ)) > max:
max = abs(expectedval(q,Id + σ))
M = expectedval(q,Id + σ)
return MLuego, se define la función que calcula las frecuencias de un modo normal para un punto del espacio recíproco dado.
def EigenfrequenciesC3(k):
"""Función que regresa las eigenfrecuencias para un vector recíproco kz dado."""
#Se define el tensor de esfuerzos de todo el sistema.
#Matriz con constante de fuerza considerando una rotación respecto al plano x - y.
φxy = pi/3
Kx2 = dot(T(φxy),dot(Kx,T(-φxy)))
#Matrices de fuerza para cada uno de los vectores en la red.
K23 = dot(U(pi),dot(Kx2,U(-pi))) #Calculada tomando 2 como origen.
K12 = dot(U(pi/3),dot(Kx2,U(-pi/3))) #Calculada tomando 1 como origen.
K13 = dot(U(-pi/3),dot(Kx2,U(pi/3))) #Calculado tomando 3 como origen.
DicMat = { 0: (K12 + K13)/m1,
1: -K12/sqrt(m1*m2),
2: -K13/sqrt(m1*m3)*exp(-1J*k*c),
3: -K12/sqrt(m1*m2),
4: (K23+K12)/m2,
5: -K23/sqrt(m2*m3),
6: -K13/sqrt(m1*m3)*exp(1J*k*c),
7: -K23/sqrt(m2*m3),
8: (K23 + K13)/m3 }
Dinteger = array([[0,1,2],
[3,4,5],
[6,7,8]])
Dm = [ [DicMat[i] for i in rw] for rw in Dinteger ]
Dynamical = asarray(np.bmat(Dm))
ω2,eigvecs = eigh(Dynamical)
return sqrt(abs(ω2))/2/pi,eigvecsSe definieron funciones que calculaba la polarización de fonón, la CCM y el pseudoescalar de momento angular dadas una frecuencia y un modo normal.
def CCM_general(ω,eigvecs):
CCM = []
for i in range(len(ω)):
eigvec = eigvecs[:,i]
q1 = sqrt(m1) * eigvec[0:3]
q2 = sqrt(m2) * eigvec[3:6]
q3 = sqrt(m3) * eigvec[6:]
Q = [q1,q2,q3]
num = 0
den = 0
for q in Q:
num += CCMMax(q)
den += vdot(q,q).real
CCM.append(1 - (num/(2*den)))
return CCM# Primero se define una función que regresa el operador de polarización de fonón.
def OperadorPseudomomento(n):
Sz = array([[0,-1J,0],
[1J,0,0],
[0,0,0]],dtype=complex)
return np.kron(np.eye(n,dtype=complex),Sz)
S = OperadorPseudomomento(3)
def Pol_general(ω,eigvecs):
Pol = []
for i in range(len(ω)):
eigvec = eigvecs[:,i]
Pol.append(expectedval(eigvec,S))
return Poldef Pz_general(ω,eigvecs):
pzlist = []
for i in range(len(ω)):
eigvec = eigvecs[:,i]
ω2 = ω[i]
q1 = sqrt(m1) * eigvec[0:3]
q2 = sqrt(m2) * eigvec[3:6]
q3 = sqrt(m3) * eigvec[6:]
Q = [q1,q2,q3]
pz = 0
for i in range(len(Q)):
pz = pz + M[i]*Q[i][2]*(X[i]*Q[i][1] - Y[i]*Q[i][0]).real
pzlist.append(pz)
return pzlistFinalmente, se muestra la gráfica de bandas para la cadena \(C_{3}\) junto con los respectivos valores de CCM en cada punto:
Kpoints = np.linspace(-pi/c,pi/c,500)
Klabels = [r"K/2",r"Γ",r"K/2"]
Kticks = concatenate([[0],[250],[500]])
BandasC3 = []
EigvecsC3 = []
CCMC3 = []
PolC3 = []
PzC3 = []
KevaluateC3 = (list( map(EigenfrequenciesC3,Kpoints) ))
for ω2,eigvec in KevaluateC3:
BandasC3.append(ω2)
EigvecsC3.append(eigvec)
for i in range(len(BandasC3)):
CCMC3.append(CCM_general(BandasC3[i],EigvecsC3[i]))
PolC3.append(Pol_general(BandasC3[i],EigvecsC3[i]))
PzC3.append(Pz_general(BandasC3[i],EigvecsC3[i]))
CCMC3 = array(CCMC3)
BandasC3 = array(BandasC3)
PolC3 = array(PolC3)
PzC3 = array(PzC3)numKC3,nbandsC3 = shape(BandasC3)
kenumC3 = arange(numKC3)
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(4.5,6))
for j in range(nbandsC3):
graficaC3 = ax.scatter(kenumC3,BandasC3.T[j],c = PolC3.T[j],cmap = "coolwarm",vmax = 1, vmin = -1,s = 5)
fig.colorbar(graficaC3)
ax.set_xticks(Kticks)
ax.set_xticklabels(Klabels,fontsize = 20);
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(4.5,6))
for j in range(nbandsC3):
graficaC3 = ax.scatter(kenumC3,BandasC3.T[j],c = CCMC3.T[j],cmap = "Reds",vmax = 1, vmin = 0,s = 5)
fig.colorbar(graficaC3)
ax.set_xticks(Kticks)
ax.set_xticklabels(Klabels,fontsize = 20);
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(4.5,6))
for j in range(nbandsC3):
graficaC3 = ax.scatter(kenumC3,BandasC3.T[j],c = PzC3.T[j],cmap = "coolwarm",vmax = 0.1, vmin = -0.1,s = 5)
fig.colorbar(graficaC3)
ax.set_xticks(Kticks)
ax.set_xticklabels(Klabels,fontsize = 20);/home/deigobau/anaconda3/envs/chiral/lib/python3.12/site-packages/matplotlib/cbook.py:1699: ComplexWarning:
Casting complex values to real discards the imaginary part
/home/deigobau/anaconda3/envs/chiral/lib/python3.12/site-packages/matplotlib/axes/_axes.py:4455: ComplexWarning:
Casting complex values to real discards the imaginary part
Ahora, hacemos el estudio de quiralidad mediante el cálculo del Pseudomomento escalar. Dado un modo \(k\), este está determinado por \[(\sum p^{z}L^{z})_{k} = -(\omega_{k})^{2}\left[\sum_{i}m_{i}c^{z}_{ki}\left(x_{i}c^{y}_{k,i} - y_{i}c^{x}_{k,i}\right)\right].\] Debajo, se define una función que calcula el valor del escalar de pseudomomento para cada valor de \(k\).
## Animación de modos normales
def Animacion(k,n,num = 4):
ωchirl,Eigenvectors = EigenfrequenciesC3(k)
ωchirl = ωchirl[n]
Eigenvectors = Eigenvectors[:,n]
A = 0.5
"""NX,NY,NZ = meshgrid( range(-num + 1,num),range(-num+1,num),range(-num+1,num))
NX = NX.flatten()
NY = NY.flatten()
NZ = NZ.flatten()
NNs = column_stack([NX,NY,NZ])"""
NZ = array( range(-num + 1,num))
NNs = column_stack([NZ])
NZ = NZ.flatten()
NNs = column_stack([NZ])
avecs1 = [avec]
xr,yr,zr = dot( NNs,avecs1 ).T
xeq1,yeq1,zeq1 = dot( NNs,avecs1 ).T + + kron(r1,ones(((2*num-1),1))).T
xeq2,yeq2,zeq2 = dot( NNs,avecs1 ).T + kron(r2,ones(((2*num-1),1))).T
xeq3,yeq3,zeq3 = dot( NNs,avecs1 ).T + kron(r3,ones(((2*num-1),1))).T
R = array([xr,yr,zr]).T
if ωchirl < 10E-3:
T = 10
else:
T = 2*pi/ωchirl # periodo
veces = 5 # número de periodos
# Posición en el tiempo cero
t0 = 0
u10 = A*exp(-1J*ωchirl*t0+1J*dot(R,array([0,0,k])))
x10 = xeq1 + Eigenvectors[0]*u10
y10 = yeq1 + Eigenvectors[1]*u10
z10 = zeq1 + Eigenvectors[2]*u10
x20 = xeq2 + Eigenvectors[3]*u10
y20 = yeq2 + Eigenvectors[4]*u10
z20 = zeq2 + Eigenvectors[5]*u10
x30 = xeq3 + Eigenvectors[6]*u10
y30 = yeq3 + Eigenvectors[7]*u10
z30 = zeq3 + Eigenvectors[8]*u10
time = linspace(t0,veces*T,endpoint=False)
Frames = []
for t in time:
u1 = A*exp(-1J*ωchirl*t+1J*dot(R,array([0,0,k])))
x1 = xeq1 + Eigenvectors[0]*u1
y1 = yeq1 + Eigenvectors[1]*u1
z1 = zeq1 + Eigenvectors[2]*u1
x2 = xeq2 + Eigenvectors[3]*u1
y2 = yeq2 + Eigenvectors[4]*u1
z2 = zeq2 + Eigenvectors[5]*u1
x3 = xeq3 + Eigenvectors[6]*u1
y3 = yeq3 + Eigenvectors[7]*u1
z3 = zeq3 + Eigenvectors[8]*u1
Frames.append( go.Frame(data=[ go.Scatter3d(x=x1.real,
y=y1.real,
z=z1.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}'),
go.Scatter3d(x=x2.real,
y=y2.real,
z=z2.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}'),
go.Scatter3d(x=x3.real,
y=y3.real,
z=z3.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}') ]) )
fig = go.Figure(
data=[ go.Scatter3d(x=x10.real,
y=y10.real,
z=z10.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}'),
go.Scatter3d(x=x20.real,
y=y20.real,
z=z20.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}'),
go.Scatter3d(x=x30.real,
y=y30.real,
z=z30.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {t:.2f}') ],
frames = Frames
)
fig.update_layout(xaxis= dict(range=[-0.5, 0.5], autorange=False),
yaxis= dict(range=[-0.5, 0.5], autorange=False),
scene={
'xaxis': {'range': [-5.2*num, 5.2*num], 'rangemode': 'tozero', 'tickmode': "linear", 'tick0': -5, 'dtick': 1},
'yaxis': {'range': [-5.2*num, 5.2*num], 'rangemode': 'tozero', 'tickmode': "linear", 'tick0': -5, 'dtick': 1},
'zaxis': {'range': [-5.2*num, 5.2*num], 'rangemode': 'tozero', 'tickmode': "linear", 'tick0': -5, 'dtick': 1},
'aspectratio': {
'x': 1,
'y': 1,
'z': 1, }},
showlegend=True,
title="Start Title",
updatemenus=[dict(type = "buttons",
buttons = [dict(label="Play",
method = "animate",
args = [None,
dict(frame = dict(duration=100,redraw=True),
transition = dict(duration=100),
fromcurrent = True,
mode = 'immediate')])])]),
fig.show()Animacion(-0.05,5,1)#Animacion(0.5,5,2)Animacion(0.05,5,1)kplot = 0.5*pi/c
nband = 0
r0 = positions.reshape((9,))
A = 2
def rt(t,k,nband):
ωchirl,Eigenvectors = EigenfrequenciesC3(k)
ωchirl = ωchirl[int(nband)]
Eigenvectors = Eigenvectors[:,nband]
rt = (r0 + A*Eigenvectors*exp(-1J*ωchirl*t)).real
r1 = rt[:3]
r2 = rt[3:6]
r3 = rt[6:]
return r1,r2,r3ωchirl,Eigenvectors = EigenfrequenciesC3(kplot)
ωchirl = ωchirl[nband]
periodo = 2*pi/ωchirl
t = linspace(0,periodo,500)
R = array(list(map(rt,t,itertools.repeat(kplot),itertools.repeat(nband))))
X1,Y1,Z1 = R[:,0].T
X2,Y2,Z2 = R[:,1].T
X3,Y3,Z3 = R[:,2].T
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,6))
ax.scatter(X1,Y1,c = t,cmap = "Blues",vmax = periodo, vmin = 0,s = 2*A)
ax.scatter(X2,Y2,c = t,cmap = "Reds",vmax = periodo, vmin = 0,s = 2*A)
ax.scatter(X3,Y3,c = t,cmap = "Greens",vmax = periodo, vmin = 0,s = 2*A)
cmap1 = matplotlib.cm.get_cmap('Blues')
cmap2 = matplotlib.cm.get_cmap('Reds')
cmap3 = matplotlib.cm.get_cmap('Greens')
rgba1 = cmap1(periodo)
rgba2 = cmap2(periodo)
rgba3 = cmap3(periodo)
ax.arrow(X1[0],Y1[0],(X1[1]-X1[0])*0.01,(Y1[1]-Y1[0])*0.01,head_width = 0.066*A,head_length = 0.1*A,fc = rgba1,ec=rgba1)
ax.arrow(X2[0],Y2[0],(X2[1]-X2[0])*0.01,(Y2[1]-Y2[0])*0.01,head_width = 0.066*A,head_length = 0.1*A,fc = rgba2,ec=rgba2)
ax.arrow(X3[0],Y3[0],(X3[1]-X3[0])*0.01,(Y3[1]-Y3[0])*0.01,head_width = 0.066*A,head_length = 0.1*A,fc = rgba3,ec=rgba3)
ax.set_ylabel("y",fontsize = 20)
ax.set_xlabel("x",fontsize = 20)/tmp/ipykernel_42974/1567948512.py:19: MatplotlibDeprecationWarning:
The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
/tmp/ipykernel_42974/1567948512.py:20: MatplotlibDeprecationWarning:
The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
/tmp/ipykernel_42974/1567948512.py:21: MatplotlibDeprecationWarning:
The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
Text(0.5, 0, 'x')
def MovimientoPlano(n):
ωchirl,Eigenvectors = EigenfrequenciesC3(0)
ωchirl = ωchirl[n]
Kpoints = np.linspace(-pi/c,pi/c,500)
# Movimiento para k distinto de cero
ωchirl,Eigenvectors = EigenfrequenciesC3(0)
ωchirl = ωchirl[n]
periodo = 2*pi/ωchirl
t = linspace(0,periodo,200)
R = array(list(map(rt,t,itertools.repeat(0),itertools.repeat(n))))
X10,Y10,Z10 = R[:,0].T
X20,Y20,Z20 = R[:,1].T
X30,Y30,Z30 = R[:,2].T
# Posición para k distinto de cero
Frames = []
for k in Kpoints:
ωchirl,Eigenvectors = EigenfrequenciesC3(k)
ωchirl = ωchirl[n]
periodo = 2*pi/ωchirl
t = linspace(0,periodo,200)
Rk = array(list(map(rt,t,itertools.repeat(k),itertools.repeat(n))))
X1,Y1,Z1 = Rk[:,0].T
X2,Y2,Z2 = Rk[:,1].T
X3,Y3,Z3 = Rk[:,2].T
Frames.append( go.Frame(data=[ go.Scatter(x=X1,
y=Y1,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {k:.2f}'),
go.Scatter(x=X2,
y=Y2,
mode="markers",
marker=dict(color="red",size=5),
name=f'time = {k:.2f}'),
go.Scatter(x=X3.real,
y=Y3.real,
mode="markers",
marker=dict(color="green",size=5),
name=f'time = {k:.2f}') ]) )
fig = go.Figure(
data=[ go.Scatter3d(x=X10.real,
y=Y10.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {k:.2f}'),
go.Scatter3d(x=X20.real,
y=Y20.real,
mode="markers",
marker=dict(color="red",size=5),
name=f'time = {k:.2f}'),
go.Scatter3d(x=X30.real,
y=Y30.real,
mode="markers",
marker=dict(color="blue",size=5),
name=f'time = {k:.2f}') ],
frames = Frames
)
fig.update_layout(xaxis= dict(range=[-2., 2.], autorange=False),
yaxis= dict(range=[-2., 2.], autorange=False),
showlegend=True,
title="Start Title",
updatemenus=[dict(type = "buttons",
buttons = [dict(label="Play",
method = "animate",
args = [None,
dict(frame = dict(duration=100,redraw=True),
transition = dict(duration=100),
fromcurrent = True,
mode = 'immediate')])])]),
fig.show()MovimientoPlano(4)